function  writeSubjData(fs)

fid = fopen(['LoudnessData.txt'],'w');
fidc = fopen(['LoudnessCorr.txt'],'w');

fprintf(fid,'Time\tSubjectiveData\tLoudness\tAweightedSPL\tLongTermLoudness\tShortTermLoudness\tFilename\tFile\tPermutation\tPreference\tFamiliarity\tParticipant\tPresOrder\tGender\tAge\tMus\n');
fprintf(fidc,'R-squaredMGBST\tLagMGBST\tP-valMGBST\tFile\nR-squaredDLM\tLagDLM\tP-valDLM\tFile\n');

tree = dataStorageTree;
l = ['ABCDEF'];
for i = 1:6
  for j = 1:6
    % [fit(i,j),maxpct(i,j), maxlag(i,j)] = writeSubjDataFile(['File' num2str(i) l(j)],fs,fid,tree);
    writeSubjDataFile(['File' num2str(i) l(j)],fs,fid,fidc,tree);
  end
end
fclose(fid);
fclose(fidc);



function [out,maxcorr, maxlag] = writeSubjDataFile(file,fs,fid,fidc,tree)

if ~exist('exportsubjdata')
  exportsubjdata = 0;
end
  
disp(file);
rownum = 1;


% Get Median Subjective Data
[tsArray,data,time]  = subjDataForFile(file, fs,  tree);
medData = median(data,2);
stdErr = std(data') / sqrt(size(data,2));

% Get SLM
dsArr = get(tree,'Name','SPL A-weighted Fast','audiofile',file);
tsSLM = getData(dsArr,1);
tsSLM = resample(tsSLM.DataObj.tsObj, tsArray{1}.DataObj.Time);

% Get ThirdOctave Spectral Centroid
dsArr = get(tree,'Name','Spectral Centroid','audiofile',file,'analyserName','CPBFFT');
tsCPBSC = getData(dsArr,1);
tsCPBSC = resample(tsCPBSC.DataObj.tsObj, tsArray{1}.DataObj.Time);

% Get ThirdOctave Spectral Centroid
dsArr = get(tree,'Name','Loudness','audiofile',file,'analyserName','LoudnessCF');
tsDLM= getData(dsArr,1);


% Get MGB Loudness
dsArr = get(tree,'Name','Long-term Loudness','audiofile',file,'analyserName','LoudnessMGB');
tsMGB1= getData(dsArr,1);
tsMGB1.DataObj.tsObj.Time = tsMGB1.DataObj.tsObj.Time/10000;



% Get MGB Loudness
dsArr = get(tree,'Name','Short-term Loudness','audiofile',file,'analyserName','LoudnessMGB');
tsMGB2= getData(dsArr,1);
tsMGB2.DataObj.tsObj.Time = tsMGB2.DataObj.tsObj.Time/10000;



% Get MGB Loudness
dsArr = get(tree,'Name','Instantaneous Loudness','audiofile',file,'analyserName','LoudnessMGB');
tsMGB3= getData(dsArr,1);
tsMGB3.DataObj.tsObj.Time = tsMGB3.DataObj.tsObj.Time/10000;


%Resample all.
tsDLM = resample(tsDLM.DataObj.tsObj, tsArray{1}.DataObj.Time);
tsMGB1 = resample(tsMGB1.DataObj.tsObj, tsArray{1}.DataObj.Time);
tsMGB2 = resample(tsMGB2.DataObj.tsObj, tsArray{1}.DataObj.Time);
tsMGB3 = resample(tsMGB3.DataObj.tsObj, tsArray{1}.DataObj.Time);


len = min([length(medData) length(tsDLM.Data) length(tsSLM.Data) length(tsCPBSC.Data)]);
dat = [tsArray{1}.DataObj.Time medData tsSLM.Data(1:len) tsCPBSC.Data(1:len) tsDLM.Data(1:len)];

% save(file, 'dat');


for j = 1:length(tsArray)
  try
    z = tsArray{j}.DataObj.DataInfo.UserData.PresOrder;
  catch
    continue;
  end
  for i = 1:len
  
    fprintf(fid,'%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%s\t%s\t%s\t%s\t%s\t%s\t%d\t%s\t%d\t%.1f\n',...
    tsArray{j}.DataObj.tsObj.Time(i),...
    ((tsArray{j}.DataObj.tsObj.Data(i)/32768)*100),...
    tsDLM.Data(i),...
    tsSLM.Data(i),...
    tsMGB1.Data(i), ...
    tsMGB2.Data(i), ...
    file, ...
    file(1:5), ...
    file(6),...
    tsArray{j}.DataObj.DataInfo.UserData.Pref,...
    tsArray{j}.DataObj.DataInfo.UserData.Fam,...
    tsArray{j}.DataObj.Name(strfind(tsArray{j}.DataObj.Name,'-')+1:end),...
    tsArray{j}.DataObj.DataInfo.UserData.PresOrder,...
    tsArray{j}.DataObj.DataInfo.UserData.Gender{1},...
    tsArray{j}.DataObj.DataInfo.UserData.Age,...
    tsArray{j}.DataObj.DataInfo.UserData.MusTraining);

  end
end
 
ts1 = medData;
ts2 = tsMGB1;
% fprintf('Length ts1 %d, Length ts2 %d',length(ts1),length(ts2))

len1 = find(isnan(medData),1,'first');
if isempty(len1), len1 = length(medData); end
len2 = find(isnan(tsDLM.Data),1,'first');
if isempty(len2), len2 = length(tsDLM.Data); end

len= min(len1, len2);
%tvec = 0:1/fs:tend;
 


% ts1 = resample(ts1,tvec);
% ts2 = resample(ts2,tvec);
% len = min(length(ts1.Data),length(ts2.Data));
% plotyy(ts1.Time(fs*2:len-fs*2),ts1.Data(fs*2:len-fs*2)/32768,ts2.Time(fs*2:len-fs*2),ts2.Data(fs*2:len-fs*2));
% pause;
% 

[val1,lag,pval]  = crosscorrelation( medData(1:len), tsMGB1.Data(1:len), 'S', 'MGB1', 20);
[m,mi] = max(val1);
fprintf(fidc,'%.3f\t%i\t%.8f\t%s\t',val1(mi), lag(mi), pval(mi), file);

[val1,lag,pval]  = crosscorrelation( medData(1:len), tsDLM.Data(1:len), 'S', 'DLM', 20);
[m,mi] = max(val1);
fprintf(fidc,'%.3f\t%i\t%.8f\t%s\n',val1(mi), lag(mi), pval(mi), file);




%% Time series analysis
% [fit,out]                   = ARXTS( ts2.Data(fs*2:len-fs*2),       ts1.Data(fs*2:len-fs*2)/32768, fs, 1,2,2);
% %plot(lag,val1,'b')
% %hold on;
% [m,mi] = max(val1);
% %plot(lag(mi),m,'ro')
% maxcorr = m;
% maxlag =lag(mi);




function  [val1,lag,pval]  = crosscorrelation(ts1, ts2, ts1name, ts2name, lags)
% [val1,lag]  = crosscorrelation(ts1, ts2, lags, ts1name, ts2aname, lags)
%
% Returns how much ts2 lags ts1, and correlation values

ts1=diff(ts1);
ts2=diff(ts2);



col = 1;
% Do SPL Pearson Correlation Analysis
% for n =lags:-1:1; 
%   [r,p,rlo,rup] = corrcoef([ts1(1:end-n) ts2(n+1:end)]); 
% %   fprintf('Lag:%d R:%.3f P:%.3f %s vs. %s \n',n,r(1,2),p(1,2), ts1name, ts2name);
%   val1(col)= r(1,2);
%   cfi(1:2,col) = [rlo(1,2) rup(1,2)];
%   lag(col)= n *-1;
%   pval(col) = p(1,2);
%   col = col + 1;
% end

% Do SPL Pearson Correlation Analysis
for n =0:lags; 
  [r,p,rlo,rup] = corrcoef([ts1(n+1:end) ts2(1:end-n)]); 
 %  fprintf('Lag:%d R:%.3f P:%.3f %s vs. %s \n',n,r(1,2),p(1,2), ts2name, ts1name);
  val1(col)= r(1,2);
  cfi(1:2,col) = [rlo(1,2) rup(1,2)];
   %plot([ts1(n+1:end) ts2(1:end-n) ])
  lag(col) = n; 
   pval(col) = p(1,2);
   col = col + 1;
%   pause;
end

[m,mi] = max(abs(val1));
fprintf('Results, %s vs. %s, Max:%.3f Lag:%d P:%.7f \n', ts1name, ts2name, val1(mi), lag(mi), pval(mi));





function [fit,out] = ARXTS(ObjectiveData, SubjectiveData, fs, lags, b, a)

% Differencing
dObjectiveData = [0; diff(ObjectiveData)];
dSubjectiveData = [0; diff(SubjectiveData)];

dat = iddata(SubjectiveData(1:end),ObjectiveData(1:end),1/fs);


m = arx(dat,[b a lags]);
[fit,out] =  compare(m,dat);
compare(m,dat)
% fprintf('Fit, %.3f\n', out);



